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Abstract 

This study investigates the unsteady flow of Powell-Eyring fluid past an inclined stretching sheet. Unsteadiness in the flow is 
due to the time-dependence of the stretching velocity and wall temperature. Mathematical analysis is performed in the 
presence of thermal radiation and non-uniform heat source/sink. The relevant boundary layer equations are reduced into 
self-similar forms by suitable transformations. The analytic solutions are constructed in a series form by homotopy analysis 
method (HAM). The convergence interval of the auxiliary parameter is obtained. Graphical results displaying the influence of 
interesting parameters are given. Numerical values of skin friction coefficient and local Nusselt number are computed and 
analyzed. 
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Introduction 

The study of boundary layer flow and heat transfer over a 
stretching sheet has gained considerable attention due to its 
numerous practical applications such as paper production, hot 
rolling, drawing of plastic films, annealing and tinning of copper 
wires and metal spinning. Wang [1] proposed the problem of 
unsteady two-dimensional boundary layer flow of liquid film on 
unsteady stretching sheet. Later Andersson et al. [2] extended 
Wang's problem for heat transfer effects by considering time- 
dependent wall temperature. Further Elbashbeshy and Bazid [3] 
investigated the thermal boundary layer in the time dependent 
flow (occupying a semi-infmite domain) over an unsteady 
stretching surface. Ishak et al. [4] studied heat transfer over an 
unsteady stretching permeable surface with prescribed wall 
temperature. Radiation effects on the flow and heat transfer over 
an unsteady stretching surface with internal heat generation were 
analyzed by Abd El-Aziz [5] . Shateyi and Motsa [6] examined the 
radiation effects on the time dependent flow of liquid film on 
unsteady stretching sheet with heat and mass transfer. They 
obtained an analytic solution of the resulting problem by 
Chebyshev pseudo-spectral collocation method. Tsai et al. [7] 
investigated the flow and heat transfer over an unsteady stretching 
surface with non-uniform heat source. Mukhopadhyay [8] 
numerically analyzed the flow over unsteady permeable stretching 
sheet with variable suction and time-dependent surface temper- 
ature. In this study, the fluid with variable viscosity and variable 
thermal conductivity was taken into consideration. Analytic 
solutions for radiation effects on mixed convection flow of Jeffrey 
fluid and heat transfer past an unsteady stretching sheet were 



provided by Hayat et al. [9] . Three dimension elastico-viscous flow 
over an unsteady stretching sheet has been discussed by Hayat 
et al. [10]. Mukhopadhyay [11] extended the work [8] for flow 
near a stagnation-point with variable free stream. MHD stagna- 
tion-point flow of an electrically conducting Casson fluid past an 
unsteady stretching surface was explored by Bhattacharyya [12]. 
Yang and Baleanu [13] investigated the fractal heat conduction 
problem. They solved by using local fractional variation iteration 
method. Yang et al. [14] presented local fractional Fourier series 
solutions for non-homogeneous heat equations arising in fractal 
heat flow with local fractional derivative. 

It has now been widely recognized that in industrial and 
engineering applications, non-Newtonian fluids are more suitable 
than Newtonian fluids. Due to the flow diversity in nature, the 
rheological features of non-Newtonian fluids cannot be captured 
by a single constitutive relationship between stress and shear rate. 
For this reason, a variety of non-Newtonian fluid models 
(exhibiting different rheological effects) are available in the 
literature [15,16]. Amongst those is the Powell-Eyring fluid [17] 
which although mathematically complex has tendency to describe 
the flow behavior at low and high shear rates. It can be used to 
formulate the flows of modern industrial materials such as 
powdered graphite and ethylene glycol. Unidirectional flow of 
Powell-Eyring fluid between parallel plates with couple stresses 
was studied by Eldabe et al. [18]. Pulsatile flow of Powell-Eyring 
fluid was examined by Zueco and Beg [19]. Homotopy 
perturbation analysis of slider bearing lubricated with Powell- 
Eyring fluid was presented by Islam et al. [20] . Three-dimensional 
flow of Powell-Eyring fluid past a wedge was discussed by Patel 
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and Timol [21]. Boundary layer flow of Powell-Eyring fluid over a 
moving flat plate was analyzed by Hayat et al. [22]. Recently 
steady flow of Powell-Eyring fluid over an exponentially stretching 
sheet was numerically investigated by Mushtaq et al. [23]. It has 
been noted that literature is scarce for unsteady flow of Powell- 
Eyring fluid. To our information, the flow and heat transfer of the 
Powell-Eyring fluid thin film over an unsteady stretching sheet are 
examined by Khader and Megahed [24]. Impact of uniform 
suction/injection in unsteady Couette flow of Powell-Eyring fluid 
is explored by Zaman et al. [25]. 

The present work considers the boundary layer flow of Powell- 
Eyring fluid over an unsteady stretching sheet. The stretching 
sheet is considered inclined. In addition the effects of radiation and 
non-uniform heat source/sink are also taken into account. 
Radiative heat transfer in the boundary layer flow is very 
important from application point of view, because the quality of 
the final product is very much dependent on the rate of heat 
transfer of the ambient fluid particles. Such radiative effects are 
also important in many non-isothermal cases whereas the heat 
generation/absorption in moving fluids is significant in the 
applications involving heat removal from nuclear fuel debris, 
underground disposal of radioactive waste material, storage of 
food stuffs, dislocating of fluids in packed bed reactors and 
several others. Similar situations prevail during the manufac- 
ture of plastic and rubber sheets where it is often necessary to 
blow a gaseous medium through the not-yet solidified material, 
and where the stretching force may be varying with time. The 
dimensionless mathematical problems are solved analytically by 
homotopy analysis method (HAM) [26—40]. Homotopy analysis 
method (HAM) is one of the most efficient methods in solving 
different type of nonlinear equations such as coupled, decou- 
pled, homogeneous and non-homogeneous. Many previous 
analytic methods have some restrictions in dealing with non- 
linear equations. For illustration, in contrast to perturbation 
method, HAM is independent of any small or large parameters 
and or the existence of auxiliary parameter provides us with a 
simple way to control and adjust the convergence region which 
is a main lack of previous techniques. Also, HAM provides us 
with great freedom to choose different initial guesses to express 
solutions of the nonlinear problem. Numerical values of 
wall velocity and temperature gradient are computed and 
examined. 

Mathematical Formulation 

We consider unsteady two-dimensional incompressible flow of 
Powell-Eyring fluid past a stretching sheet. The sheet makes an 
angle a with the vertical direction. The x - and y-axes are taken 
along and perpendicular to the sheet respectively. In addition the 
effects of thermal radiation and non-uniform heat source/sink are 
considered (see Fig. 1). The Cauchy stress tensor in Powell-Eyring 
fluid is given by [17]: 
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In the above expressions t is the time, v = ( \i j p ) is the 
kinematic viscosity, k is the thermal conductivity of the fluid, p is 
the fluid density, T is the fluid temperature, c p is the specific heat, 
g 0 is the acceleration due to gravity, fi T is the volumetric 

coefficient of thermal exponential, q r = 00 -r— [36-38] is 

jk Oy 

the linearized radiative heat flux, k* is the mean absorption 
coefficient, cr* is the Stefan-Boltzmann constant, q" is the non- 
uniform heat generated ( q" > 0) or absorbed ( q" < 0) per 
unit volume. The non-uniform heat source/ sink, q" is modeled by 
the following expression [39-40]. 



q „ = ^M [A{Ts _ Tx)f+iT _ T ^ 



(4) 



in which A and B are the coefficient of space and temperature- 
dependent heat source/sink, respectively. Here two cases arise. 
For internal heat generation A>0 and B>0 and for internal heat 
absorption, we have A<0 and B<0. 

bx 

The surface velocity is denoted by u s ( x, t) = — 

whereas the surface temperature T s ( x, t) = + 
bx 2 

T re f^2 — ^ ^ _ at)~^l 2 . Here b (stretching rate) and a are 

1 . Also T re f is a 



positive constants having dimension time 
constant reference temperature. We note that the temperature of 
stretching sheet is larger than the free stream temperature T CJ0 . 
The boundary conditions are taken as follows: 



u = u s (x, t),v = 0, T = T s (x, t) at y = 0, 



t/— ►(), 7"— >• T f as y^co. 



(5) 



dui 1 j 
T i7 = A* ^ r- ^ sinn 



1 dui 

Cdxj 



where (I is the viscosity coefficient, (3 and G are the material fluid 
parameters. The boundary layer equations comprising the balance 
laws of mass, linear momentum and energy can be written as [19— 
26] 



Introducing 
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Figure 1. Physical model and coordinate system. 

doi:1 0.1 371 /journal. pone.01 0321 4.g001 
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Figure 2. The ft -curves for the velocity field. 

doi:1 0.1 371 /journal. pone.01 0321 4.g002 
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-1.5 -1.0 -0.5 0.0 

Figure 3. The // -curves for the temperature field. 

doi:1 0.1 371 /journal. pone.01 0321 4.g003 



Eq. (1) is identically satisfied and Eqs. (2)— (5) become 

/ = 0,/'= 1,0=1 at*/ = 0, 

( 1 + T)f"'-f' 2 +ff" - Tpf" 2 r -z(f'+\ nf") + GO cos cc = 0, (7) />-<), 0-0 as tj- oo, (9) 

4 \ / j \ where prime denotes differentiation with respect to rj , f is the 

1-| — R) 0"+ PvlfO'- 2f'0 e(36 + rjO') ) + Af + B6 = 0, dimensionless stream function, 0 is the dimensionless temperature 

/ , Q , and the dimensionless numbers are 



3 J V 2 



Table 1 . Convergence of series solutions for different order of approximations when a = n/4, (3 = 0.5, r = 0.2, R = 0.2, i e = 0.6, G = 0.3, 
Pr=1.0, A = B = 0.1, "r>= -0.8 and f\= -0.7. 





Order of approximation 


- /"( 0) 


- 6 '( 0) 


1 


1.03515 


1.33250 


5 


1 .04402 


1.35252 


10 


1 .04401 


1.35252 


15 


1 .04401 


1.35252 


20 


1 .04401 


1.35252 


30 


1 .04401 


1.35252 



doi:1 0.1 371/journal.pone.01 03214.t001 
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Figure 4. Influence of oc on the velocity field. 

doi:1 0.1 371 /journal. pone.01 0321 4.g004 
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Figure 5. Influence of G on the velocity field. 

doi:1 0.1 371 /journal. pone.01 0321 4.g005 
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Figure 6. Influence of p on the velocity field. 

doi:1 0.1 371 /journal. pone.01 0321 4.g006 



Local Nusselt number Nu x is defined as 




2 4 6 8 

Figure 7. Influence of e on the velocity field. 

doi:1 0.1 371 /journal. pone.01 0321 4.g007 
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Figure 8. Influence of r on the velocity field. 

doi:1 0.1 371 /journal.pone.01 0321 4.g008 



where Re x = is the local Reynolds number. 

v 

Solution Methodology 

Most of the problems occurring in the field of science and 
engineering are non-linear. Specifically most of the problems 
encountered in fluid mechanics are highly non-linear. To find the 
exact solution of these non-linear problems is very difficult and 
some times even impossible. Thus several numerical and analytical 
techniques have been developed to solve such kind of problems. 



Among these HAM is the most used analytical technique. 
Convergent series solutions of non-linear equations are obtained. 

Homotopy analysis method 

HAM was proposed by means of homotopy, a fundamental 
concept of topology. Two functions are said to be homotopic if one 
function can be deformed continuously into the other function. If 
f\ and f2 are two continuous maps from the topological space X 
into the topological space Y then f\ is homotopic to fa if there exist 
a continuous map F 




Figure 9. Influence of R on the velocity field. 

doi:1 0.1 371 /journal.pone.01 0321 4.g009 
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£ = 0.5 : a = - 4. G = 03, T = 0.2 : e = 0.6 : Pr = 1.0. A = B = 0.1 




R = 0.0.R = 0.3.R = 0.6.R = 0.9 



— i L_ 



-I L 



—i L 



i: 
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Figure 10. Influence of R on the temperature field. 

doi:1 0.1 371 /journal. pone.01 0321 4.g01 0 



F : Ix[0,l]^7 



such that for each x e X 



F(x,G)=f x {x),F(x,\)=f 2 (x) 



It should be noted that there is a great freedom to choose initial 
guess and auxiliary linear operator Beside such a great freedom 
there are some fundamental rules which direct us to choose the 
mentions parameters in more efficient way. Therefore, initial 
guesses for the velocity and temperature fields are taken in such a 
way that they satisfy the boundary conditions given in Eq. (9). And 
we choose linear operator specified in Eq. (13) that must satisfy the 
properties given in Eq. (14). 



The map F is called homotopy between fa and fa. 



(12) 



£ = 0.5 : a = nf4, G = 0.3. T = 0.2 : € = 0.6 : R = O.2.. A = B = 0.1 




Pr=1.0.Pr=1.3.Pr=1.6.Pr=1.9 



Figure 11. Influence of Pr on the temperature field. 

doi:1 0.1 371 /journal.pone.01 0321 4.g01 1 
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? = 0.5 : <r = N/4, G = 0.3 : Pr = 1.0 : 6 = 0.<5 : R = 0.2 : A = B = 0.1 




r»o.o,r«o.5,r«i.5,r-2.o 



Figure 12. Influence of F on the temperature field. 

doi:1 0.1 371 /journal.pone.01 0321 4.g01 2 



£f = f -f, £o = e -9, 

subject to the properties 

£y(d + C 2 e" + Cjc-') =0, £ e (C 4 e" + Cjc"') =0, (13) 

where C\ (i= 1-5) are the constants. 



The deformation problems subjected to zeroth order 

(l-^^Pj-ZoWl^^P)]. (14) 

(1 -P)£e[0(t,;p) - floW] =/>M%5/>)L (15) 



I = 0.5, <* = a/4, T = 0.2 : Pr = 1.0. € = 0.6. R = O.2.. A = B = 0.1 




G = 0.0.G = 0.5.G=1.5.G = 2.5 



Figure 13. Influence of G on the temperature field. 

doi:1 0.1 371 /journal.pone.01 0321 4.g01 3 
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? = 0.5 : a = - 4 : r = 0.2 : Pt = 1.0. € = 0.6 : R = 0.2 : G = 0.3 : B = 0.1 




A. = -0.3, A = -0.1. A = 0.2, A = 0.3 



Figure 14. Influence of A on the temperature field. 

doi:1 0.1 371 /journal.pone.01 0321 4.g01 4 



/(0; p) = 0,/(0; p) = 1 ,/(oo ; p) = 0, 0(0; p) = 1 , 0(oo ; p) = 0. (16) 



/(»?; 1) =/(»?), %; 1) = %). 



If p e [0,1] indicates the embedding parameter, hf and kg the 
non-zero auxiliary parameters then the nonlinear differential 
operators N/ and are given by 



It is noticed that when p varies from 0 to 1 then f(fj',p) and 
6 (r\ ; p) approach from fo( rj ) , 6 o( v\ ) to f(f]) and 
6 (rj ) . The series of / and 0 through Taylor's expansion are 
chosen convergent for p = 1 and thus 



N/[fti/;/>)] = (l+r: 



5 2 /0/;/0 fa/fa; /0 



d/7 2 



P ^ 5>? 2 J £ ^ 5>j 2^ <3>? 2 

+ G0(>?; />) cos a, 



/0/)=/ofa)+ £/»0/),/m(»/)= 



1 gTM (20) 



m=Uri)+ E WU>?)= ^^%=o- (21) 



N 9 [%; p),f(m P)Y- 



£ 39(t]-p) + t] 



se{n;p) 

drj 



+A 



CI] 

drj 



CI] 



(18) 



+B9(r,;p). 



The resulting problems at order can be presented in the 
following forms 



£f[fm(ri)-x m f m -i(n)]=fifRf m (ri), 



£e[G m {n)-x m e m -x{n)\=fifRe m {n), 



(22) 



(23) 



We have for p = 0 and p = 1 the following equations 

AmO)=Mn)MmO)=0o(fi), (19) 



/«(0) =/ M (0) =/ m (oo) = 0,„(O) = 0 m (oo) 
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Figure 15. Influence of B on the temperature field. 

doi:1 0.1 371 /journal.pone.01 0321 4.g01 5 




Figure 16. Influence of e on the temperature field. 

doi:1 0.1 371 /journal.pone.01 0321 4.g01 6 
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Table 2. Comparison between numerical solution Tsai et. al. [7] and HAM solution in a special case when 

a = p = s = r = G= R= 0. 




Pr B 


A 


Present study Tsai et. [7] 


1.0 - 1.0 


0.0 


- 1.71094 - 1.710937 


- 2.0 


- 1.0 


- 2.36788 


2.0 - 1.0 


0.0 


- 2.25987 


- 2.0 


- 1.0 


- 2.486000 - 2.485997 


doi:10.1371/journal.pone.0103214.t002 






x m = 0, m< 1, 




0m(n) = Q* m (n) + C A e* + C 5 e-\ (27) 



x m = 1, m>0. 



The general solutions (f m , 0 m ) comprising the special 
solutions (f m * , 0 * ) are 



(26) 



Convergence of the homotopy solutions 

It is now a well established argument that the convergence of 
series solutions (22) and (23) depends upon the auxiliary 
parameters h . The admissible range of values of hf and ho (for 
some fixed values of parameters) lie along the line segment 



parallel to hf and ho 



axes. For example in Figs. 2 and 3 



the permissible range of values of hf and ho are 



Table 3. Values of heat transfer characteristics at wal 


- 


0 '( 0) for different emerging parameters when t> = 


-0.8 and ^ =-0.7. 




a r p g 


G 


R Pr A 




0.0 






1.35702 


7U/6 






1.35798 


TU/3 






1 .34926 


tt/4 0.0 






1 .69881 


0.4 






1.72556 


0.7 






1.74109 


0.9 






1 .74984 


0.0 






1.71555 


0.5 






1.71319 


0.9 






1.71114 


0.0 






1.10162 


0.4 






1.61674 


0.6 






1.71319 




0.0 




1 .69868 




0.5 




1.72227 




0.8 




1.73515 






0.0 


1 .54046 






0.3 


1.71319 






0.6 


2.00303 






1.2 


1.91058 






1.5 


2017721 






1.9 


2049323 






-0.1 


1 .80783 






0.0 


1 .76059 






0.1 


1.71319 
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- 1.26 < h f < - 0.24 and - 1.4 < h 0 < - 0.25 respectively 
when 8 = 0.6. This series solutions converge for the whole region 
of rj when hf= —0.9 and hg = —0.8. Table 1. shows the 
convergence of HAM solution for different order of approxima- 
tions. It is clear from this table that 10 th order of approximations 
are sufficient for convergent solutions up to six decimal place. 

Results and Discussion 

This section examines the effects of different physical param- 
eters on the velocity and temperature fields. Hence Figs. (4, 5, 6, 7, 
8, 9, 10, 11, 12, 13, 14, 15, 16) are plotted. Fig. 4 elucidates the 
behavior of inclination angle a on the velocity and the boundary 
layer thickness, a = 0 shows the corresponding velocity profiles in 
the case of a vertical sheet for which the fluid experiences the 
maximum gravitational force. On the other hand when a changes 
from 0 to 71 / 2 i.e. when the sheet moves from vertical to 
horizontal direction, the strength of buoyancy force decreases and 
consequently the velocity and the boundary layer thickness 
decrease. Fig. 5 indicates that velocity field /' is an increasing 
function of G. This is because a larger value of G accompanies a 
stronger buoyancy force which leads to an increase in the l x — 
component of velocity. The boundary layer thickness also 
increases with an increase in G. Variation in f with an increase 
in /? can be seen from Fig. 6. It is noticed that /' decreases and 
boundary layer thins when ft is increased. Influence of unsteady 
parameter e on the velocity field is displayed in Fig. 7. Increasing 
values of £ indicates smaller stretching rate in the x - direction 
which eventually decreases the boundary layer thickness. Inter- 
estingly the velocity increases by increasing £ at sufficiently large 
distance from the sheet. Variation in the x-component of velocity 
with an increase in the fluid parameter T can be described from 
Fig. 8. In accordance with Mushtaq et al. [25], the velocity field 
/' increases with an increase in T . 

Radiation effects on the velocity and temperature distributions 
are perceived from Figs. 9 and 10. An increase in R enhances the 
heat flux from the sheet which gives rise to the fluid's velocity and 
temperature. Wall slope of the temperature function therefore 
increases with an increase in R. Fig. 1 1 portrays the effect of 
Prandtl number on the thermal boundary layer. From the 
definition of Pr given in Eq. (10), it is obvious that increasing 
values of Pr decreases conduction and enhances pure convection 
or the transfer of heat through unit area. That is why temperature 
and the thermal boundary layer thickness decrease with an 
increase in Pr . This reduction in the thermal boundary layer 
accompanies a larger heat transfer rate from the sheet. Temper- 
ature profiles for different values of T are shown in Fig. 12. It is 
seen that temperature 6 is an increasing function of T . Fig. 1 3 
indicates that an increase in the strength of buoyancy force due to 
temperature gradient decreases the temperature and the thermal 
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boundary layer thickness. Influence of heat source/ sink parame- 
ters on the thermal boundary layer are presented in Figs. 14 and 
15. As expected the larger heat source (corresponding to A>0 and 
B>0) rises the fluid's temperature above the sheet. While the non- 
uniform heat sink corresponding to A<0 and B<0 can contribute 
in quenching the heat from stretching sheet effectively. Fig. 16 
depicts that temperature 6 is a decreasing function of the 
unsteady parameter £ . 

Table 2 shows comparison of present work with Tsai et al. in a 
special case. A very good agreement is found between the results of 
wall temperature gradient. Table 3 shows the effect of embedded 
parameters on heat transfer characteristics at the wall — 6 ' ( 0) . 
Since in the present case the sheet is hotter than the fluid i.e 
T w > Too thus heat flows from the sheet to the fluid and hence 
6 ' ( 0) is negative. From this table we observe that with an 
increase in a , /? and R the wall heat transfer rate | 0 ' ( 0) | 
decreases. However it increases when Y , £ and Pr are 
increased. 

Conclusions 

This article addressed the radiation effects in the unsteady 
boundary layer flow of Powell-Eyring fluid past an unsteady 
inclined stretching sheet with non-uniform heat source/sink. 
Homotopy analysis method (HAM) was used to obtain approx- 
imate analytic solutions of the governing nonlinear problem. The 
important findings are listed below. 

1) The strength of gravitational force can be varied by changing 
the inclination angle a which the sheet makes with the 
vertical direction. The velocity decreases with an increase in 
a . 

2) Velocity field /' and temperature 6 are decreasing function 
of the unsteady parameter £ . 

3) Velocity increases and temperature decreases when the fluid 
parameter T is increased. 

4) Increase in the radiation parameter R enhances the heat flux 
from the plate which increases the fluid's velocity and 
temperature. 

5) The analysis for the case of viscous fluid can be obtained by 
choosing r = ft = 0. Further the results for horizontal 
stretching sheet are achieved for a = n /2. 
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